GENERALIZED GALERKIN VARIATIONAL INTEGRATORS 



MELVIN LEOK 

Abstract. We introduce generalized Galerkin variational integrators, which are a natural generalization 
of discrete variational mechanics, whereby the discrete action, as opposed to the discrete Lagrangian, is 
the fundamental object. This is achieved by approximating the action integral with appropriate choices 
of a finite-dimensional function space that approximate sections of the configuration bundle and numeri- 
cal quadrature to approximate the integral. We discuss how this general framework allows us to recover 
higher-order Galerkin variational integrators, asynchronous variational integrators, and symplectic-energy- 
momentum integrators. In addition, we will consider function spaces that are not parameterized by field 
values evaluated at nodal points, which allows the construction of Lie group, multiscale, and pseudospectral 
variational integrators. The construction of pseudospectral variational integrators is illustrated by applying 
it to the (linear) Schrodinger equation. G-invariant discrete Lagrangians are constructed in the context of 
Lie group methods through the use of natural charts and interpolation at the level of the Lie algebra. The 
reduction of these G-invariant Lagrangians yield a higher-order analogue of discrete Euler-Poincare reduc- 
tion. By considering nonlinear approximation spaces, spatio-temporally adaptive variational integrators can 
be introduced as well. 
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1. Introduction 

We will review some of the previous work on discrete mechanics and their multisymplectic generalizations 
(see, for example, Marsden et al. [1998, 2001]), before introducing a general formulation of discrete mechanics 
that recovers higher-order variational integrators (see, for example, Marsden and West [2001]), asynchronous 
variational integrators (see, for example, Lew et al. [2003]), as well symplectic-energy- momentum integrators 
(see, for example, Kane et al. [1999]). 

While discrete variational integrators exhibit desirable properties such as symplecticity, momentum preser- 
vation, and good energy behavior, it does not address other important issues in numerical analysis, such 
as adaptivity and approximability. Generalized variational integrators are introduced with a view towards 
addressing such issues in the context of discrete variational mechanics. 
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By formulating the construction of a generalized variational integrator in terms of the choice of a finite- 
dimensional function space and a numerical quadrature scheme, we are able to draw upon the extensive 
literature on approximation theory and numerical quadrature to construct variational schemes that are 
appropriate for a larger class of problems. Within this framework, we will introduce multiscale, spatio- 
temporally adaptive, Lie group, and pseudospectral variational integrators. 

1.1. Standard Formulation of Discrete Mechanics. The standard formulation of discrete variational 
mechanics (see, for example, Marsden and West [2001]) is to consider the discrete Hamilton's principle, 

5§ d = 0, 

where the discrete action sum, : Q n+1 — > M, is given by 

rt-l 

§ d (q ,qi, q n ) = L d {q l , q l+1 ). 

i=0 

The discrete Lagrangian, Ld : Q x Q —> R, is a generating function of the symplectic flow, and is an 
approximation to the exact discrete Lagrangian, 

LT C \qo,qi) = / L(q 01 (t),q 01 (t))dt, 
Jo 

where goi(0) = qo, qoi{h) = qi, and goi satisfies the Euler-Lagrange equation in the time interval (0, h). The 
exact discrete Lagrangian is related to the Jacobi solution of the Hamilton-Jacobi equation. The discrete 
variational principle then yields the discrete Eulei — Lagrange (DEL) equation, 

D 2 Ld{qa, qi) + D 1 L d (q 1 ,q 2 ) = 0, 

which yields an implicit update map (<Zoj<Zi) l— > (ftj'fe) that is valid for initial conditions (qo,qi) that as 
sufficiently close to the diagonal of Q x Q. 

1.2. Multisymplectic Geometry. The generalization of the variational principle to the setting of partial 
differential equations involves a multisymplectic formulation (see, for example, Marsden et al. [1998, 2001]). 
Here, the base space X consists of the independent variables, which are denoted by (x°, . . . , x n ), where x° 
is time, and x 1 , . . . ,x n are space variables. 

The independent or field variables, denoted (q 1 , . . . , q m ), form a fiber over each space-time basepoint. The 
set of independent variables, together with the field variables over them, form a fiber bundle, tt : Y — > X, 
referred to as the configuration bundle. The configuration of the system is specified by giving the field 
values at each space-time point. More precisely, this can be represented as a section of Y over X, which is a 
continuous map q : X — > Y, such that tt o q = l x . This means that for every x G X, q(x) is in the fiber over 
x, which is tt^ 1 (x). 

In the case of ordinary differential equations, the Lagrangian is dependent on the position variable, 
and its time derivative, and the action integral is obtained by integrating the Lagrangian in time. In the 
multisymplectic case, the Lagrangian density is dependent on the field variables, and the derivatives of the 
field variables with respect to the space-time variables, and the action integral is obtained by integrating the 
Lagrangian density over a region of space-time. 

The analogue of the tangent bundle TQ in the multisymplectic setting is referred to as the first jet 
bundle J X Y ', which consists of the configuration bundle Y, together with the first derivatives of the field 
variables with respect to independent variables. We denote these as, 

dxJ 

for i = 1, . . . , m, and j = 0, . . . , n. 
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We can think of J X Y as a fiber bundle over X. Given a section q : X — > Y, we obtain its first jet 
extension, j 1 q : X — ► J X Y, that is given by 

j\{x°, . . . , X n ) = (x°, . . . , X n , q\ . . . , q m , q\ u . . . , q m J , 

which is a section of the fiber bundle J X Y over X. 

The Lagrangian density is a map L : J X Y — > il n+1 (X), and the action integral is given by 



S(q) - f Ltfq), 
J X 



and then Hamilton's principle states that 

(55 = 0. 

We will see in the next subsection how this allows us to construct multisymplectic variational integrators. 

f.3. Multisymplectic Variational Integrator. We introduce a multisymplectic variational integrator 
through the use of a simple but illustrative example. We consider a tensor product discretization of (1 + 1)- 
space-time, given by 



ij+i- 















At 

Ax 





















Xi+l 





Tj-1- 

Xi-1 

and tensor product shape functions given by 

1 

<Pij(x,t) = 

Xi Xi-\-\ 

We construct the discrete Lagrangian as follows, 

Ld(qi,j,qi+i,j,q%,j+i,qi+i,j+i) = / ^f^fE ■ E! • 1a,b^Pa,b)) , 

J[x i ,x i+1 ]J[t J ,t J+1 ] V \^a =l ^b=j )) 

where qij ~ q(iAx, jAt). 

Consider a variation that varies the value of qij , and leaves the other degrees of freedom fixed, then we 
obtain the following discrete Euler-Lagrange equation, 

DiL c i(q i j,q i+ ij,q i j + i,q i+ ij + i)+D2L c i(q i _ij,q i j,q i _ij + i,q i j + i) 

+DsL c i(q i j_i,q i+ ij_i,q i j,q i+ ij) 

+DiL<i{qi-\ y i-i,qi y i-i,qi-\ y i,qi,i) = 0, 
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In general, when we take variations in a degree of freedom, we will have a discrete Euler-Lagrange equation 
that involves all terms in the discrete action sum that are associated with regions in space-time that overlap 
with the support of the shape function associated with that degree of freedom. 



There are a few essential observations that go into constructing a general framework that encompasses 
the prior work on variational integrators, asynchronous variational integrators, and symplectic-energy- 
momcntum integrators, while yielding generalizations that allow the construction of multiscale, spatio- 
temporally adaptive, Lie group, and pseudospectral variational integrators. 

The first is that a generalized Galerkin variational integrator involves the choice of a finite-dimensional 
function space that discretizes a section of the configuration bundle, and the second is that we approximate 
the action integral though a numerical quadrature scheme to yield a discrete action sum. The discrete 
variational equations we obtain from this discrete variational principle are simply the Karush-Kuhn-Tucker 
(KKT) conditions (see, for example, Nocedal and Wright [1999]) with respect to the degrees of freedom that 
generate the finite-dimensional function space. 

To recap, the choices which are made in discretizing a variational problem are: 

(1) A finite-dimensional function space to represent sections of the configuration bundle. 

(2) A numerical quadrature scheme to evaluate the action integral. 

Given these two choices, we obtain an expression for the discrete action in terms of the degrees of freedom, 
and the KKT conditions for the discrete action to be stationary with respect to variations in the degrees of 
freedom yield the generalized discrete Euler-Lagrange equations. 

Current variational integrators are based on piecewise polynomial interpolation, with function spaces that 
are parameterized by the value of the field variables at nodal points, and a set of internal points. By relaxing 
the condition that the interpolation is piecewise, we will be able to consider pseudospectral discretizations, 
and by relaxing the condition that the parameterization is in terms of field values, we will be able to 
consider Lie group variational integrators. By considering shape functions motivated by multiscale finite 
elements (see, for example, Hou and Wu [1999], Efendiev et al. [2000], Chen and Hou [2003]), we will obtain 
multiscale variational integrators. And by generalizing the approach used in symplectic-energy-momentum 
integrators, and considering nonlinear approximation spaces (see, for example, DeVore [1998]), we will be 
able to introduce spatio-temporally adaptive variational integrators. 

As we will see, there is nothing canonical about the form of the discrete Euler-Lagrange equations, or 
the notion that the discrete Lagrangian is a map Ld ■ Q x Q — > K. These expressions arise because of the 
finite-dimensional function space which has been chosen. 

2.1. Special Cases of Generalized Galerkin Variational Integrators. In this subsection, we will show 
how higher-order Galerkin variational integrators, multisymplectic variational integrators, and symplectic- 
energy-momentum integrators are all special cases of generalized Galerkin variational integrators. 

Higher-Order Galerkin Variational Integrators. In the case of higher-order Galerkin variational in- 
tegrators, we have chosen a piecewise interpolation for each time interval [0, ft], that is parameterized by 
control points q$, . . .q^, corresponding to the value of the curve at a set of control times = c?o < d\ < 
. . . < ds-i < d s = 1. The interpolation within each interval [0, h] is given by the unique degree s polynomial 
q d (t; q%, h), such that q d {d v h) = g%, for v = 0, . . . , s. 

By an appropriate choice of quadrature scheme, we can break up the action integral into pieces, which we 
denote by 



2. Generalized Galerkin Variational Integrators 
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dih d 2 h d s - 2 h d s -ih h t 

Figure 1: Polynomial interpolation used in higher-order Galerkin variational integrators. 



If we further require that the piecewise defined curve is continuous at the node points, we obtain the following 
augmented discrete action, 



v=0,. 



N-l 



N-2 



The discrete action is stationary when 



i=0 



dqt 



i=0 



^ 9 +i + ~ 



as 

dqj 



We can identify the pieces of the discrete action with the discrete Lagrangian, L d : Q x Q — > R, by setting 

(2.1) M«i,9i+i) = si(«n. 

where, = qf = qi+i, and the other qf's arc defined implicitly by the system of equations 



(2.2) 



dq: 



for v = 1, . . . , s — 1. Once we have made this identification, we have that 



gi+l 

-£>iM%+i>%+2) = -^4— W+i) 

= A, 

= D 2 L d (q l ,q l+ i), 
from which we recover the discrete Euler-Lagrange equation, 

DiL d (q i+1 ,qi +2 ) + D 2 L d (qi 7 q i+1 ) = 0. 
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The DEL equations, together with the definition of the discrete Lagrangian, given in Equations 2.1 and 2.2, 
yield a higher-order Galerkin variational integrator. As we have shown, it is only because the interpolation 
was piecewise that we were able to decompose the equations into a DEL equation, and a set of equations 
that define the discrete Lagrangian in terms of conditions on the internal control points. 

Multisymplectic Variational Integrators. Recall the example of a multisymplectic variational integrator 
we introduced in §1.3, where we used tensor product linear shape functions with localized supports to 
discretize the configuration bundle. The discrete action is then given by 

N N 

§d ( tia,b}a=0 M-l) = I / L (j 1 ( / 'InJrrn.l, 

\ b=0,...,N 



= / i (j 1 (EE ? »-^ 1 

J[x ,xm] J [Wiv] a= o b=0 

M-1N-1 r r N N 

EE/ / iWEEw. 

»=0 j=0 J {*i,Xi+i] • y [*j>*j+i] a=06=0 



3- 

M-1N-1 i+1 j+l 

= EE/ / ^MEE^ 

where we first decomposed the integral into pieces, and then used the local support of the shape functions 
to simplify the inner sums over q a ,b^fa,b- As before, it is due to the local support of the shape functions that 
we can express the discrete action as 

M-l JV-l 

§><i = E E Ld(qi,j,qi+i,j,qi,j+i,qi+i,j+i), 

i=0 j=0 



where 



/ 


f L ( 


/ [Xi,X i + 1 ] ■ 


'[tj,t j+1 ] 



Ld{qi,j,qi+i,j,qi,j+i,qi+i,j+i) = / / L [j 1 [ E a=i E 6= -- 9a,b<p a ,b 



This localized support is also the reason why the discrete Euler-Lagrange equation in this case consists 
of four terms, since to each degree of freedom, there are four other degrees of freedom that have shape 
functions with overlapping support. In the case of ordinary differential equations, this was two. In general, 
for tensor product meshes of (n + l)-space-times, with tent function shape functions, the number of terms 
in the discrete Euler-Lagrange equation will be 2" +1 . In contrast, for pseudospectral variational integrators 
with m spatial degrees of freedom per time level, and k degrees of freedom per piecewise polynomial in time, 
each of the m(k — 1) discrete Euler-Lagrange equations will involve m(k — 1) terms. 

This is simply a reflection of the fact that shape functions with compact support yield schemes with 
banded matrix structure, whereas pseudospectral and spectral methods tend to yield fuller matrices. The 
payoff in using pseudospectral and spectral methods for problems with smooth or analytic solutions is due 
to the approximation theoretic property that these solutions are approximated at an exponential rate of 
accuracy by spectral expansions. 

Symplectic-Energy-Momentum Integrators. In the case of symplectic-energy-momentum integrators, 
the degrees of freedom involve both the base variables and the field variables. We will first derive a second- 
order symplectic-energy momentum integrator, and in the next section, we will derive a higher-order general- 
ization. We choose a piecewise linear interpolation for our configuration bundle. Each piece is parameterized 
by the endpoint values of the field variable qf, q}, and the endpoint times h®, h}, and we approximate the 
action integral using the midpoint rule. To ensure continuity, we require that, qj = qf + i, and hj = hf +1 . 

Remark 2.1. The approach of allowing each piecewise defined curve to float around freely , and imposing the 
continuity conditions using Lagrange multipliers was used in Lall and West [2003] to unify the formulation 



GENERALIZED GALERKIN VARIATIONAL INTEGRATORS 



7 



of discrete variational mechanics and optimal control. Through the use of a primal-dual formalism, discrete 
analogues of Hamiltonian mechanics and the Hamilton- Jacobi equation were also introduced. 



This yields the following discrete action, 



JV-l 



i=0 



q" + ql q}-Q°i 



2 'h\-h\ 



JV-2 



N-2 



i=0 



i=0 



To simplify the expressions, we define 



hi = hi - hi 
qf + oj 



Then, the variational equations are given by 

0=L [q iH ,q i+i )-h~(q iH , qiH ) ±-q i+{ 



L (q i+ i,q^ 



dL 



d-q\ q ^ q ^)^ q ^ 



=hi 
=hi 



dL 
dq~ 



l \ 1 dL ( \ 1 

<9£ / . \ 1 dL ( . \ 1 



— A^, 



h 1 -h° 

If we define the discrete Lagrangian to be 

Ld(qi, qi+i,hi) = hiL 

and the discrete energy to be 

Ed(Qi,Qi+i, hi) = -D 3 L d (qi,q i+1 , hi) 

,'qi + q i+ i q i+ i - q, 

and identify points as follows, 
we obtain 



for i = 0, . . . , N - 2, 



for i = 1,...,N- 1, 
for i = 0, . . . , N - 2, 

for i = l,...,iV- 1, 

for i = 0, . . . , N - 2, 
for i = 0, . . . , N - 2. 



hi 



Qi = 



hi 



dL fqi + q i+1 q l+ i - q, t \ 1 q i+1 - q { 



94 



= 9°, 



- hU = hi 



=E d (qi,q i+1 ,hi) - u u 
= - q i+1 , hi) + u>i-i, 

=D 2 L d (qi, qi+i, hi) - A i; 
=DiL d {qi, qt+i, hi) + Aj_i, 
After eliminating the Lagrange multipliers, we obtain the conservation of discrete energy equation, 

Ed(qi,qi+i, hi) = E d (q i+ i,q i+ 2, hi+i), 



for i 


= 0,. 


.,7V- 


2- 


for i 


= 1,- 


.,7V- 


1, 


for i 


= 0,. 


.,7V- 


2, 


for i 


= 1,. 


.,7V- 


1. 
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and the discrete Euler-Lagrangc equation, 

D 2 L d (q i ,q l+ i,h l ) + DiL d (q i+1 , q i+2 , h i+ i) = 0. 

This recovers the results obtained in Kane et al. [1999], but the derivation is new. In §5, we will see how 
this derivation can be generalized to yield a higher-order scheme. 

Discrete Action as the Fundamental Object. The message of this section is that the discrete action is 
the fundamental object in discrete mechanics, as opposed to the discrete Lagrangian. In instances whereby 
the shape function associated with individual degrees of freedom have localized supports, it is possible 
to decompose the discrete action into terms that can be identified with discrete Lagrangians. While this 
approach might seem artificial at first, we will find that in discussing pseudospectral variational integrators, 
it does not make sense to break up the discrete action into individual pieces. 

3. Lie Group Variational Integrators 

In this section, we will introduce higher-order Lie group variational integrators. The basic idea behind all 
Lie group techniques is to express the update map of the numerical scheme in terms of the exponential map, 

9i = 3oexp(£ i) , 

and thereby reduce the problem to finding an appropriate Lie algebra element £oi £ 0, such that the update 
scheme has the desired order of accuracy. This is a desirable reduction, as the Lie algebra is a vector space, 
and as such the interpolation of elements can be easily defined. In our construction, the interpolatory method 
we use on the Lie group relies on interpolation at the level of the Lie algebra. 

For a more in depth review of Lie group methods, please refer to Iserles et al. [2000]. In the case 
of variational Lie group methods, we will express the variational problem in terms of finding Lie algebra 
elements, such that the discrete action is stationary. 

As we will consider the reduction of these higher-order Lie group integrators in the next section, we will 
chose a construction that yields a G-invariant discrete Lagrangian whenever the continuous Lagrangian is 
G-invariant. This is achieved through the use of G-equivariant interpolatory functions, and in particular, 
natural charts on G. 

3.1. Galerkin Variational Integrators. We first recall the construction of higher-order Galerkin varia- 
tional integrators, as originally described in Marsden and West [2001]. Given a Lie group G, the associated 
state space is given by the tangent bundle TG. In addition, the dynamics on G is described by a La- 
grangian, L : TG — > R. Given a time interval [0, h], the path space is defined to be 

C(G) = C([0, h],G) = {g : [0, h] -► G | g is a C 2 curve}, 

and the action map, & : C(G) — > M, is given by 

6(g) = I'* L(g(t),g(t))dt. 
Jo 

We approximate the action map, by numerical quadrature, to yield & s : C([0, h], G) — > M, 

s 

& s {g) = hY J hL{g{c l h),g{c t h)), 

i=l 

where q € [0, 1], i = 1, . . . , s are the quadrature points, and bj are the quadrature weights. 
Recall that the discrete Lagrangian should be an approximation of the form 

LJg ,g\,h) « ext 6(g) . 

1 ' geC([o,h],G)MO)=9oMh)=9i W 
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If we restrict the extremization procedure to the subspace spanned by the interpolatory function that is 
parameterized by s + 1 internal points, ip : G s+1 — > C([0, h], G), we obtain the following discrete Lagrangian, 

L d (g ,9i)= ext G(T<p(g"; •)) 

g iy eG;g u =go;g s =gi 

s 

= ext h^2b 1 L(Tip(g u ;c l h)). 

g"£G;g°=g ;g 3 =gi T - f 

The interpolatory function is G-equivariant if 

<p(gg u ;t) = gv{g v ;t). 

Lemma 3.1. If the interpolatory function ip(g u ;t) is G-equivariant, and the Lagrangian, L : TG — ► K., is 
G-invariant, then the discrete Lagrangian, Ld : G x G — > M, gwen &y 



Ld(go,gi)= ext h'S^b,L(TLp(g v ;c t h)), 



g" £G;g°=g ;g s =gi 



is G-invariant. 
Proof. 

s 

Ld(ggo>ggi) = ext h V^L(T^(^;c;/i)), 

g u £G;g°=ggo\g s =ggi 

2—1 

5 

ext hY / b l L(T^(gg--c l h)), 

g"£g L G\g"=g ;g s =gi — ' 
s 

= ext /i> biL(TL g ■ Tip(g l/ ;cih)), 

g l, eG;g°=g ;g 3 =g 1 

i=l 

s 

= ext h bjL(Tip(g v ; Cjh)), 

g u eG;g =g ;g B =gi T - : 
i=l 

= L d (g ,gi), 

where we used the G-equivariance of the interpolatory function in the third equality, and the G-invariancc 
of the Lagrangian in the forth equality. □ 

Remark 3.1. While G-equivariant interpolatory functions provide a computationally efficient method of 
constructing G-invariant discrete Lagrangians, we can construct a G-invariant discrete Lagrangian (when 
G is compact) by averaging an arbitrary discrete Lagrangian. In particular, given a discrete Lagrangian 
Ld ■ Q x Q — > R, the averaged discrete Lagrangian, given by 



L d {q ,qi) = 7^77 / L d (gq a ,gqi)dg 

M Jg£G 



is G-equivariant. Therefore, in the case of compact symmetry groups, a G-invariant discrete Lagrangian 
always exists. 

3.2. Natural Charts. Following the construction in Marsden et al. [1999], we use the group exponential 
map at the identity, exp e : g — + G, to construct a G-cquivariant interpolatory function, and a higher-order 
discrete Lagrangian. As shown in Lemma 3.1, this construction yields a G-invariant discrete Lagrangian if 
the Lagrangian itself is G-invariant. 
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In a finitc-dimcnsional Lie group G, cxp e is a local diffcomorphism, and thus there is an open neighborhood 
U C G of e such that expj 1 : U — > u C g. When the group acts on the left, we obtain a chart tjj g : L g U — > u 
at g <E G by 

Lemma 3.2. The interpolatory function given by 

is G-equivariant. 
Proof. 

ip(gg v ;rh) = ^ g \o ) { y Y^ v=0 ^9 B (gg V %A T )) 

= L gg o cx Pe fel^eH(99T 1 (99 v ))lAT)) 

= Lg^g* ( E„ =0 ex P^ oL (g r 1 (9 v )LA t )) 

= M> 1 (£l =0 ^^)M r )) 

= L g ip{g v ;Th). □ 

Remark 3.2. In the proof that ip is G-equivariant, it was important that the base point for the chart should 
transform in the same way as the internal points g v . As such, the interpolatory function will be G-equivariant 
for a chart that it based at any one of the internal points g v that parameterize the function, but will not be 
G-equivariant if the chart is based at a fixed g G G. Without loss of generality, we will consider the case 
when the chart is based at the first point go . 

We will now consider a discrete Lagrangian based on the use of interpolation in a natural chart, which is 
given by 

s 

L d (g , gi )= ext ftYJ ^L(T^({^}« =0 ; ah) . 

g l, eG;g°=ga:,g 3 =g a g± j=1 

To further simplify the expression, we will express the extremal in terms of the Lie algebra elements 
associated with the i^-th control point. This relation is given by 

and the interpolated curve in the algebra is given by 

s 

^ V ;rh) = J2ClAr), 

K=0 

which is related to the curve in the group, 

9(9 v ;Th) = g exp(£(ip go (g v );Th)). 

The velocity £ = g~ Y g is given by 

i(rh) = g^girh) = ±J2cIAt)- 

Using the standard formula for the derivative of the exponential, 

T s exp = T e L eMi) ■ dcxp adi; , 
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where 

oo 

dexp„, = 

n=0 

we obtain the following expression for discrete Lagrangian 



(n + iy: 



Ld(go,gi)= ext h y^b l L[L go exp(^(c l h)), 

5-e fl ;C°=0;C»= 1 /' S o(9i) ~[ V 



T, 



\^p{i{ Ci h))L go ■ T e L 0X p( ?(Ci/l )) • dexp ad ^^ h) (£(cih))J . 

More explicitly, we can compute the conditions on the Lie algebra elements for the expression above to be 
extremal. This implies that 

s 

Ld(9o,9i) = h^2biL^L go exp(£(c;/i)), T exp{i{c . h)) L go • T e L exp{( . {c . h)) ■ dexp^^ (£(c;/i))) 



i=i 

with £° = 0, £ s = ipg (gi), and the other Lie algebra elements implicitly defined by 
' dL 

— (c l h)T cxp{({c . h)) L go ■ T e L cxp{({c . h)) ■ dexp ad5(Cj)j) l v , s (a) 
1 dL 



Q = hJ2b l 



i=l 



+ ^■^■( c ^) T o X p(?(c l / l )) i o X p(c( c ,ft)) • T e L e xp(|(c 4 h)) • ddex Pad 4(C4h) l v . s { Cl 



for v = 1, . . . ,s — 1, and where 

oo 



n=0 v ' 

This expression for the higher-order discrete Lagrangian, together with the discrete Euler-Lagrange equation, 

DiLd{go,gi) + D 1 L d (g ll g 2 ) = 0, 
yields a higher-order Lie group variational integrator. 

4. Higher-Order Discrete Euler-Poincare Equations 

In this section, we will apply discrete Euler-Poincare reduction (see, for example, Marsden et al. [1999]) 
to the Lie group variational integrator we derived previously, to construct a higher-order generalization of 
discrete Euler-Poincare reduction. 

4.1. Reduced Discrete Lagrangian. We first proceed by computing an expression for the reduced discrete 
Lagrangian in the case when the Lagrangian is G- invariant. Recall that our discrete Lagrangian uses G- 
equivariant interpolation, which, when combined with the G-invariancc of the Lagrangian, implies that the 
discrete Lagrangian is G-invariant as well. We compute the reduced discrete Lagrangian, 

Idig^gi) = L d (g ,gi) 

= Ldie.g^gx) 



s 

ext h/^ biL[L e exp(£(cih)), 

;? o =0 . e=log(s -i Sl) ^ V 

T eX p(£( Ci h))Le ■ T e L cxp(s(Cth)) • dexp ad?(cjh) (£(cj/i))) 

S 

ext h V]Wd exp(t(ah)),T e L eMi{c . h)) • dexp ad ac . h) (i(c t h))) . 
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Setting £° = 0, and £ s = log(g gi), we can solve the stationarity conditions for the other Lie algebra 
elements {£"}*=i using the following implicit system of equations, 



= h^h 



dL 

{cih)T e L eMi{c . h)) ■ dexp ad?(cth) l Vta (ci) 



dg 



1 dL ~ 

-f l ^{cih)T^L eMi{c . h)) • ddexp ad?(cs(i) L, s {ci) 



where v = 1 , . . . , s — 1. 

This expression for the reduced discrete Lagrangian is not fully satisfactory however, since it involves 
the Lagrangian, as opposed to the reduced Lagrangian. If we revisit the expression for the reduced discrete 
Lagrangian, 

s 

ld(go 1 9i)= ext ftV6 i i(exp($(c i /i)),r e L exp( f (c4h)) -dexp^ (^(cift))) , 

we find that by G-invariance of the Lagrangian, each of the terms in the summation, 

L^exp(^(cih)),T e L exp{s{c . h)) ■ dexp ad£(cih) (£(c;/i))) , 

can be replaced by 

^( dcx Pad £(Cih) (£M))) . 

where I : q — > K is the reduced Lagrangian given by 

l{n) = L(L g -ig,TL g -ig) = L(e,n), 

where rj = TL g -\g G g. 

From this observation, we have an expression for the reduced discrete Lagrangian in terms of the reduced 
Lagrangian, 

s 

ld(9o 1 9i)= ext /iVM(dex Pad £(c . M (£M)))- 

As before, we set £° = 0, and £ s = log(<7 _1 <7i), and solve the stationarity conditions for the other Lie algebra 
elements using the following implicit system of equations, 



= hJ2 b * 



dl ~ 

— (c t h) ddex Pa(W) l V)S ( Ci ) 



where v = 1, . . . ,s — 1 . 



4.2. Discrete Euler Poincare Equations. As shown above, we have constructed a higher-order reduced 
discrete Lagrangian that depends on 

./fcfc+i = 9kg k+i- 

We will now recall the derivation of the discrete Euler-Poincare equations, introduced in Marsden et al. 
[1999]. The variations in fkk+i induced by variations in g k , g k +i are computed as follows, 

Sfkk+i = -g^&gkgk-tgk+i +9k ls 9k+i 

= TR hk+i(-9k ls 9k + Ad /fcfc+1 g k+1 Sg k+1 ) . 
Then, the variation in the discrete action sum is given by 

N-l 

<5§> = ^ l'd(fkk+l)S fkk+l 
k=0 



GENERALIZED GALERKIN VARIATIONAL INTEGRATORS 13 
N-l 



= E l 'd(fkk+i)TR fkk+1 (-g k 1 5g k + Ad fkk+1 g k+1 Sg k+1 ) 

k=0 
N-l 

= E [Uh-ik)TR fk _ lk kd h _ lk -l' d (f kk+1 )TR fkk+1 ] $ k , 



fe=i 



with variations of the form = 5g k - In computing the variation of the discrete action sum, we have 
collected terms involving the same variations, and used the fact that t? = #jv = 0. This yields the discrete 
Euler-Poincare equation, 

l' d (fk-ik)TR fk _ lk Ad fk _ lk -l' d (f kk+1 )TR fkk+1 = 0, k = 1, . . • , N - 1. 

For ease of reference, we will recall the expressions from the previous subsection that define the higher-order 
reduced discrete Lagrangian, 



i=i 

where 



ld(fkk+l) = /lE^( deX Pad, (Cih) (CM))) 
i=l 

s 



K = 

and 

e° = o, 

r = l0g(/fcfc+l) , 

and the remaining Lie algebra elements are defined implicitly by 



dl ~ 
— ( Cl /i)ddexp ad5(ct(i) l v>a (a) 



= hJ2 b 

i=l 

for i/ = 1 , . . . , s — 1, and where 

oo n 

ddex P- = E (^)T • 

n=0 v ; 

When the discrete Euler-Poincare equation is used in conjunction with the higher-order reduced discrete 
Lagrangian, we obtain the higher-order Eulet — Poincare equations. 



■5. Higher-Order Symplectic-Energy-Momentum Variational Integrators 

In this section, we will generalize our new derivation of the symplectic-energy-momentum preserving 
variational integrators (see, Kane et al. [1999]) to yield integrators with higher-order accuracy. 

As before, we consider a piecewise interpolation, with both the control points in the field variables, q", 
and the endpoints of the interval, h®,hj, as degrees of freedom. The continuity conditions for this function 
space are q* — q® +1 , and h\ = h? i+1 . Then, we have that the discrete action is given by 

N-l s 

= ^( h l- h i)T, b ^i(c j (hl-h^q^Mc j (hl-h^,q^) 

i=0 j=l 

JV-2 JV-2 

-E w-«?+i)-E w *(^-^+i)> 

i=0 i=0 
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where 



qi (T(hl-h°); q ») = J2>liks(T) 



1=0 



Ur{h\ - /*?);#) = -^-^ £ qfl K , s {r). 



K = 



To simplify the expressions, we define hi = h\ — hf. Then, the variational equations are given by 



s s dL 1 

= ^2b j L(q i (c j h i ),q i (c j h i )) - hi^bj — icjh^—q^Cjhi) -un 

3 = 1 3 = 1 ° q 1 

s s dL 1 

= -^b j L(q l (c :j h i ),q i (cjh i )) + hi^bj — (cjhi)—qi(cjhi)+GJi-i, 



3 = 1 

S 

3 = 1 

s 

3 = 1 

s 

3 = 1 

Qi =Qi+i, 
h 1 -h° 

n i — n i+li 

We can eliminate the Lagrange multipliers, to yield 

=^2b J L(q l (c : jh l ),q l (c : jh t )) -^bj — ^jh^q^Cjh 



3 = 1 

dL ~ 1 dL ~ 

-Q^(cjh t )l s , s (c 3 ) + — — {c hi)l StS {c ) 

dL , , N r , \ 1 dL . ~ . ' 

-Q-(Cjhi)lo,s( c j) + J^-Qq-^ 11 *) 10 '^) 

9L , , . - , . 1 dL ~ . 



hi 



3 = 1 



3 = 1 



'y]bjL(q i - 1 (cjhi- 1 ),qi- 1 (cjhi-i)) 



3 = 1 



dL 



y^bj — (cj/ti_i)gi_i(cj/tj_i), 

3 = 1 



3 = 1 



9L , . .r , , 1 dL ~ 

— (Cjhi)l StS ( C] ) + ■j^-g^( C 3 h i) l sA C j) 



+ hi-! b 3 

3 = 1 



dL - 1 dL ~ ' 

— (cjhi-ijlo^Cj) + — - — {cjhi-^lo^Cj) 



% S =9?+i, 
ft 1 -7>° 



-g^icjhi^sicj) + — — {cjhjl^icj) 



for i = 0, . . . , N - 
for i = 1,...,N- 
for i = 0, . . . , N - 

for i = 1,...,N- 

for i = 0, . . . , N - 
v = l,...,s- 

for i = 0, . . . , N - 
for i = 0, . . . , N - 



for i = 1,...,N- 1 



for i = l,...,iV- 1 

for i = 0, . . . , N - 1 
v = 1, . . . , s — 1 

for i = 0, . . . , N - 2 
for i = 0, . . . , N - 2 
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If we define the discrete Lagrangian as follows, 

s 

L d (q l ,q l+1 ,h i ) = h t ^ bjL(qi(cjhi), qi{cjhi)), 
3=1 

where 

s 

q l {Th i -q v t ) = 5^«ff K , a (r), 

K=0 
1 S - 

Qi(rhi;q^) = — ^ <?f Z~ k , s (t), 
and = q^, q? = q\, and the remaining terms were defined implicitly by 

s 

i=i 

then the equations reduce to the following, 

d 

E d (qi,q l+1 ,hi) = ~ — [L d (q t ,q l+1 ,h l )} 7 

Ed(qi,Qi+l,hi) = E d (q i+ i,q i+ 2,h i+ i), 

= D 2 L d (qi, q i+1 , hi) + DiL d (q i+1 ,q i+2 , h i+1 ). 
which is a higher-order symplectic- energy -momentum variational integrator. 

Solvability of the Energy Equation. It should be noted that the discrete energy conservation equation is 
not necessarily solvable, in general, particularly near stationary points. This issue is discussed in Kane et al. 
[1999], Lew et al. [2004], and can be addressed by reformulating the discrete energy conservation equation 
as an optimization problem that chooses the time step by minimizing the discrete energy error squared. 
Clearly, reformulating the discrete energy conservation equation yields the desired behavior whenever the 
discrete energy conservation equation can be solved, while allowing the computation to proceed when discrete 
energy conservation cannot be achieved, albeit with a slight energy error in that case. This does not degrade 
performance significantly, since instances in which discrete energy conservation cannot be achieved are rare. 

6. Spatio-Temporally Adaptive Variational Integrators 

As is the case with all inner approximation techniques in numerical analysis, the quality of the numerical 
solution we obtain is dependent on the rate at which the sequence of finite-dimensional function spaces 
approximates the actual solution as the number of degrees of freedom is increased. 

For problems that exhibit shocks, nonlinear approximation spaces (see, for example, DeVore [1998]), as 
opposed to linear approximation spaces, are clearly preferable. Adaptive techniques have been developed in 
the context of finite elements under the name of r-adaptivity and moving finite elements (see, for example, 
Baines [1995]), and has been developed in a variational context for elasticity in Thoutireddy and Ortiz [2003]. 
The standard motivation in discrete mechanics to introduce function spaces that have degrees of freedom 
associated with the base space is to achieve energy or momentum conservation, as discussed in §5, or Kane 
ct al. [1999], Oliver et al. [2004]. However, if the solution to be approximated exhibits shocks, nonlinear 
approximation techniques achieve better results for a given number of degrees of freedom. 

In this section, we will sketch the use of regularizing transformations of the base space, to achieve a 
computational representation of sections of the configuration bundle that will yield more accurate numerical 
results. 



dL , , ,~ , . 1 dL . ~ . 
-r^{c hi)l VtS (cj) + — — {cjhi^sicj) 
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Consider the situation when we are representing a characteristic function using piecewise spline interpo- 
lation. We show in Figure 2, the difference between linear and nonlinear approximation of the characteristic 
function. 



J 



(a) Using equispaced nodes (b) Using adaptive nodes 



Figure 2: Linear and nonlinear approximation of a characteristic function. 



When the derivatives of the solution vary substantially in a spatially distributed manner, we obtain 
additional accuracy, for a fixed representation cost, if we allow nodal points to cluster near regions of high 
curvature. It is therefore desirable to consider variational integrators based on function spaces that are 
parameterized by both the position of the nodal points on the base space, as well as the field values over the 
nodes. 

This is represented by having a regular grid for the computational domain 1Z, which is then mapped to 
the physical base space X, as shown in Figure 3. 



\ 






X / / 


\ \ 



Figure 3: Mapping of the base space from the computational to the physical domain. 
The sections of the configuration bundle factor as follows, 



Y 



K 




The mapping tp : 1Z — ► X results in a regularized computational representation q : 1Z — > Y of the original 
section q : X — > Y. The relationship between the discrete section of the configuration bundle and its 
computational representation is illustrated in Figure 4. 
The action integral is then given by 

S(q) = [ Ltfq) = [ L{fq)\T>y\. 
Jx Jtz 
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Thus, even though q : X — > Y may exhibit shocks, the computational representation we work with, q : 1Z — > 
Y, is substantially more regular, and consequently, a numerical quadrature scheme in 1Z applied to 

/ L[j l <t)\T>ip\ 

is significantly more accurate than the corresponding numerical quadrature scheme in X applied to 

/ Ltfq). 

Jx 

As such, spatio-temporally adaptive variational integrators achieve increased accuracy by allowing the ac- 
curate representation of shock solutions using an adapted free knot representation, while using a smooth 
computational representation to compute the action integral. 

7. Multiscale Variational Integrators 

In the work on multiscale finite elements (MsFEM), introduced and developed in Hou and Wu [1999], 
Efendiev et al. [2000], Chen and Hou [2003], shape functions that are solutions of the fast dynamics in the 
absence of slow forces are constructed to yield finite element schemes that achieve convergence rates that 
are independent of the ratio of fast to slow scales. 

In constructing a multiscale variational integrator, we need to choose finite-dimensional function spaces 
that do a good job of approximating the fast dynamics of the problem, when the slow variables are frozen. 
In addition, we require an appropriate choice of numerical quadrature scheme to be able to evaluate the 
action, which involves integrating a highly-oscillatory Lagrangian. In this section, we will discuss how to go 
about making such choices of function spaces and quadrature methods. 

We will start with a discussion of the multiscale finite element method, to illustrate the importance of 
a good choice of shape functions in computing solutions to problems with multiple scales. After that, we 
will walk through the construction of a multiscale variational integrator for the case of a planar pendulum 
with a stiff spring. Finally, we will discuss how we might proceed if we do not possess knowledge of which 
variables, or forces, are fast or slow. 
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7.1. Multiscale Shape Functions. We will illustrate the idea of constructing shape functions that are 
solutions of the fast dynamics by introducing a model multiscale second-order elliptic partial differential 
equation given by 

V-a(x/e)Vu e (x)=f(x), 

with homogeneous boundary conditions. In the one-dimensional case, we can solve for the solution analyti- 
cally, and it has the form 



u £ (x) 



F{y) 



dy 



r 1 F (v) d y 

JO a(y/e) a y 
Jo a 



dy 



a{yhY 



where F(x) — f(y)dy. If we have nodal points at {xj}^ , then the appropriate multiscale shape functions 
to adopt in this example is to use shape functions that are solutions of the homogeneous problem at the 
element level. These shape functions ipj satisfy 

im ( a ( x l e )-ik l f'i) = °> for x *-i < x < x i+^ 
\<Pi{xi-i) = 0; <Pi{xi+i) = 0; ip\{xi) = 1. 

And they have the explicit form given by 



¥>i(*) 




ds 


-1 


r 


ds 


a(s/e)_ 




a(s/e) 


ds 


-1 


" rx i+1 
Jx 


ds 


a(s/e)_ 




a(s/e) 



x G , Xi\ , 

x G (xi : Xi+i] , 
otherwise . 



It can be shown that this will yield a numerical scheme that solves exactly for the solution at the nodal 
points. 

As an example, we will compute the analytical solution for a(x) = 1+0 QSshi^Tra) 1 f( x ) = a;2 ' an< ^ e = 0-025. 
This is illustrated in Figure 5(a). What is particularly interesting is to compare the zoomed plot of the 
exact solution and the multiscale shape function over the same interval, shown in Figures 5(b) and 5(c), 
respectively. The multiscale finite element method is able to achieve excellent results because the multiscale 
shape functions are able to capture the fast dynamics well. In the next subsection, we will discuss how this 
insight is relevant in the construction of multiscale variational integrators. 

7.2. Multiscale Variational Integrator for the Planar Pendulum with a Stiff Spring. As was 

shown previously, a shape function that captures the fast dynamics of a multiscale problem is able to achieve 
superior accuracy when used for computation. While this idea has primarily been used for problems with 
multiple spatial scales, it is natural to consider its application to a problem with multiple temporal scales, 
such as the problem of the planar pendulum with a stiff spring, as illustrated in Figure 6. 

We will use this example to illustrate the issues that arise in constructing a multiscale variational inte- 
grator. The variables are q = (a, 9), where a is the spring extension, and 8 is the angle from the vertical. 
The Lagrangian is given by 

L{a, 0, a, 9) = — {x + y ) - mgy - -a , 



where 



x = (I + a) sin0, 

y = — (I + a) cos 9, 

x = a sin 9 + 8(1 + a) cos 9, 

y = —a cos 9 + 9(1 + a) sin 9. 
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Figure 5: Comparison of the multiscale shape function and the exact solution for the elliptic problem. 
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Figure 6: Planar pendulum with a stiff spring 



The Hamilton's equations for the planar pendulum with a stiff spring are 

Pa 



a = 
m 

A Pe 



m(/ + a) 2 ' 
p a = — ka + gmcosd + m(l + a)9 , 
pe = —gm(l + a) sin6>. 

The timescale arising from the mass-spring system is 2-rry/m/k. The timescale arising from the planar 



pendulum system is 2n\/T/g. The ratio of timescales is given by e = y/mg/kl. 

Multiscale Shape Function. In this problem, the fast scale is associated with the stiff spring, and if we 
set the slow variable 9 = 0, we obtain the equation 

.. Pa k 

a = — = a, 

m m 

which has solutions of the form 

a(t) = a sin ^y/k/mtj + a\ cos ^y/k/mtj . 

We will now consider a well-resolved simulation of this system using the odel5s stiff solver from Matlab, 
with parameters m = 1, g = 9.81, k = 10000, I — 1, giving a scale separation of e — 0.0313. The simulation 
results are show in Figure 7. 

Clearly, if we wish to choose time steps that do not resolve the fast oscillations in a, but do resolve the slow 
oscillations in 9, it would be desirable to include sin ^yjk/mtj , and cos (^/k/mtj in the finite-dimensional 
function space used to interpolate a(t). 

Evaluating the Discrete Lagrangian. Since we have chosen time steps that do not resolve the fast 
oscillations, it follows that over the interval [0,h], the Lagrangian will oscillate rapidly as well. In com- 
puting the discrete Lagrangian, it is therefore necessary to ensure that this highly-oscillatory integral is 
well-approximated. 

It is conventional wisdom, in numerical analysis, that the numerical quadrature of highly-oscillatory 
integrals is a challenging problem requiring the use of many function evaluations. Recently, there has been 
a series of papers, Iserles [2003a, b, 2004], Iserles and N0rsett [2004], that provides an analysis of Filon-type 
quadrature schemes that provide an efficient and accurate method of evaluating such integrals. The method 
is applicable to weighted integrals as well, but we will summarize the results from §3 of Iserles [2003a] 
restricted to unweighted integrals, and refer the reader to the original reference for an in-depth discussion 
and analysis. 
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(a) a(t) 



(b) 9(t) 





(c) a(t) (zoomed) 



(d) 0.5 cos (Vfc/mt) 



Figure 7: Comparison of the multiscale shape function and the exact solution for the planar pendulum with 
a stiff spring. 



The Filon-type method aims to evaluate an integral of the form 

h[f] = f h f(x)e^dx = h [ H f{hx)e^dx. 
Jo Jo 

Given a set of distinct quadrature points, ci < c 2 < • • • < c v in [0, 1], the Filon-type quadrature method is 
given by 



Q^[f} = hJ2h^hLu)f( Cl h), 



i=i 



where 



h(ihuj) = / U(x)e in " x dx, 
Jo 

and U are the Lcgendrc polynomials. Here, we draw attention to the fact that the quadrature weights are 
dependent on hu>. 

If the quadrature points correspond to Gauss-Christoffel quadrature of order p, then the error for the 
Filon-type method is given as follows. 

0(h p+1 ), if huj < 1 
0{h v ), if huj = <D{l) 

0{h V+1 /{hLu)), if huj > 1 
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0(h v+1 /{hij) 2 ), if huj > 1, Cl = 0, C = 1. 

Clearly, for highly-oscillatory functions, the last case, which corresponds to the Lobatto quadrature points, is 
most desirable. Thus, it is appropriate to use the Filon-Lobatto method to evaluate the discrete Lagrangian 
in our case. 

Discrete Variational Equations. As we discussed previously, instead of looking for stationary solutions 
to the discrete Hamilton's principle in polynomial spaces, we will consider solutions that are piecewise of the 
form 

q(t;{pj},LU,a ,ai) = o Pj '^) (! + a o sm M) + a i cos ( w *)) • 

This function space approximates the highly-oscillatory nature of the solution well, in contrast to a 
polynomial function space, thereby avoiding approximation-theoretic errors. The degrees of freedom in 
this function space are {pj} 1 w, a , and a\. Since there are no distinguished degrees of freedom that are 
responsible for the endpoint values of the curve, we need to impose continuity at the nodes using a Lagrange 
multiplier. The augmented discrete action is given by 

N-l r h 

§ rf = E / {p}},^>oXi))^ 

N-2 

+ A i (gi(/i;{pj},w%4,ai)-g i+1 (0;{pj +1 },^+ 1 ,4 +1 ,ai +1 )), 

i=0 

where each of the integrals are evaluated using the Filon-Lobatto method. Taking variations with respect 
to the degrees of freedom yields an update map, 

which gives the multiscale variational integrator for the planar pendulum with a stiff spring. 

7.3. Computational Aspects. Multiscale variational integrators have the advantage of directly accounting 
for the contribution of the fast dynamics, thereby allowing the scheme to use significantly larger time-steps, 
while maintaining accuracy and stability. It is possible to take advantage of knowledge about which of the 
variables, or forces, are fast or slow, by using a low degree polynomial and oscillatory functions for the fast 
variables, and a higher-order polynomial for the slow variables. In the absence of such information, it is 
appropriate to use a function space with both polynomials and oscillatory functions, and apply it to all the 
variables. 

Recall that the Filon-type method has quadrature coefficients that depend on the frequency. As such, the 
initial fast frequency has to be estimated numerically using a fully resolved computation for a short period 
of time. Since both the function space and the quadrature weights depend on the fast frequency u, the 
resulting scheme is implicit and fairly nonlinear, and as such, it may be expensive for large systems. 

8. PSEUDOSPECTRAL VARIATIONAL INTEGRATORS 

The use of spectral expansions of the solution in space are particularly appropriate for highly accurate 
simulations of the evolution of smooth solutions, such as those arising from quantum mechanics. We will 
introduce pseudospectral variational integrators, and consider the Schrodingcr equation as an example. 

In particular we will adopt the tensor product of a spectral expansion in space, and a polynomial expansion 
in time. For example, we could have an interpolatory function of the form 

N/2 

#r,(r + 0At) = -L ]T' e <fc *((l-T)fii+Tt>i +1 ), 

k=-N/2 



GENERALIZED GALERKIN VARIATIONAL INTEGRATORS 



23 



which is the tensor product of a discrete Fourier expansion in space, and linear interpolation in time. Here, 
the Yl notation denotes a weighted sum where the terms with indices ±7V/2 are weighted by 1/2, and the 
other terms are weighted by 1. See page 19 of Trefethen [2000] for a discussion of why this is necessary to 
fix an issue with derivatives of the interpolant. 

The degrees of freedom are given by v k , which are the discrete Fourier coefficients. We will later see 
how such an interpolation can be applied to the Schrodinger equation. The action integral can be exactly 
evaluated for this class of shape functions, as we will see below. 

It is straightforward to generalize the pseudospectral approach we present in this section to a spectral 
variational integrator, with discrete Fourier expansions in space for periodic domains, or Chebyshev expan- 
sions in space for non-periodic domains, and Chebyshev expansions in time. This will however result in all 
the degrees of freedom on the space-time mesh being coupled, and is therefore substantially more expensive 
computationally than the pseudospectral method. The payoff for adopting the spectral approach is spectral 
accuracy, which is accuracy beyond all orders. 



8.1. Variational Derivation of the Schrodinger Equation. Let H be a complex Hilbert space, for 
example, the space of complex- valued functions tp on R 3 with the Hermitian inner product, 



where the overbar denotes complex conjugation. We will present a Lagrangian derivation of the Schrodinger 
equation, following worked example 9.1 on pages 568-569 of Jose and Saletan [1998]. 
Consider the Lagrangian density C given by 



where H : H — > H is given by 
which yields 



Hip = -—V 2 iP + ViP, 
2m 



if i 



Pi 2 



c{ytp) = —{ipip - ipip} - — Vtp • W> - Vtptp. 



2 1 2m 

We take tp, tp as independent variables, and compute, 



l J Ldt 'J 
-I 



( dC ^ dC ^ dC m -\ fdC tl dC tl dC _ , N 
\dtp Qijj dV?/> / dtp oVi/j 



fdC_ d_dC _ 
\dtp dtg^p 




dC 



dC 



dtp dt dip dVtp 



dtp 



d 3 xdt 



d A xdt 



tp ) 5ip ■ 



jip-ViP + -iP + —V 2 ib)6iP 



d A xdt, 



where we integrated by parts, and neglected boundary terms as the variations vanish at the boundary of 
the space-time region. Since the variations are arbitrary, we obtain the nonrelativistic (linear) Schrodinger 
equation as a result, 



mtp = {-^ + v}tp 



We note that the Lagrangian density is invariant under the internal phase shift given by 

ip i-> e le tp, tp i-» e~ le tp. 
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The space part of the multi-momentum map is given by 



dC dC — ih 2 — — 



and the time part is given by 



dip dip 2 



The norm of the wavefunction is automatically preserved by variational integrators, since the norm is a 
quadratic invariant. 

8.2. Pseudospectral Variational Integrator for the Schrodinger Equation. Consider a periodic do- 
main [0, 2tt], discretized with a discrete Fourier series expansion in space, and a linear interpolation in time. 
Let N be an even integer, then, our computation is done on the following mesh, 

71" 2tt 
I . . . . . . . . . . 

Xi X 2 X N / 2 X N -1 X N 

This implies that the grid spacing is given by 

h= N- 

The interpolation is given by 

N/2 

iP(x,(r + l)At) = ±- Y! e^((l-r)vi+Tv l + 1 ), 

k=-N/2 
N/2 

V,(*,(T + Z)Ai) = -^ e ik *(v l k +1 -v l k ), 

k=-N/2 



k=-N/2 



N/2 
1 v^' 



ip(x,(T+i)A t )=— y: 



ikx + l _ f.l s 



2irAt 

k=-N/2 

and the discrete Fourier transformation is given by 

JV 



e [v, -v 



k > 



2ir 



for k — —N/2 + 1, . . . , N/2, and V-n/2 = Vn/2- Recall that 
where H : TL — > H is given by 

h 2 

Hip = V 2 ip + Vip. 

2m 

Furthermore, the potential V is expressed using a discrete Fourier expansion, 



N/2 



2ir 

k=-N/2 
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In addition, we will need to introduce a normalization condition, so as to eliminate trivial solutions of the 
partial differential equation. The normalization condition is 



J ° \ k=-N/2 ) \ k=-N/2 ) k=-N/2 



which is enforced using a Lagrange multiplier. Discrete Action for the Schrodinger Equation. The 

discrete action in the space-time region [0, 2ir] x [I At, (I + I) At] is given by 



(/+l)At r 2-i7 



I At JO 

(i+l)At /-27T 



Ctf^dxdt + hil-^u^)) 



Jim Jo 



JO 



ih 



ih ■ - - h 2 
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2 2 771 



JV/2 



N/2 



dxdt + \i 1 1 - — E 44 
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N/2 
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k=-N/2 



N/2 



- \h E' ^(a-^* 



k=-N/2 

N/2 

■At ^ 



-ikx 13+1 ~Z ' 



k=-N/2 



2irAt 



k=-N/2 



At dxdr 



1 /-27T 1-2 



' (-*")«- ((l-rrt + r<») 

u \ k=-N/2 , 

( N/2 \ 

y k=-N/2 ) 

(N/2 \ / N/2 \ 

^ E' ^ i E' ^(d-^+^l 

k=-jV/2 } \ m=-N/2 } 



1 /-27r 

Jo 



JV/2 

^ E' e-*»((l-r)fi 

n=-N/2 
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i r n/2 

k=-N/2 



dr 



~ E" ((1 - r)v{ + r^)((l - + r^ 1 ) 

fe=-JV/2 



1 

2^ 



-l 



JV/2+n 



E' E' (K-™((l-r)^+re)((l-*+rii +1 ) 

n=-N/2 m=-N/2 
N/2 N/2 

+ E' E' ( W(l - r)v l rn + rv l +'){{\ - r)¥ n + r^ 1 )) 

n=0 m=n-N/2 
\ k=-N/2 ) 



N/2 

ifl TT^'I 



,,,,, W/2 

ft 2 /c 2 Ai ^» 



k=-N/2 
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E E 
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+ A M E" 
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where we used the fact that 



Jo 



o • 



for k E Z, and we define as a weighted sum where the terms with indices ±N/2 are weighted by 1/2, 
and as a weighted sum where the terms with indices ±N/2 are weighted by 1/4. We should note that 
using the same approach, it would be possible to exactly evaluate the action integral for the class of tensor 
product shape functions with a discrete Fourier expansion in space, and a polynomial expansion in time. In 
particular, a similar approach is valid in exactly evaluating the action integral when we use shape functions 
that are spectral in both space and time. 

Discrete Euler Lagrange Equations. We are now in a position to compute the discrete Euler-Lagrangc 
equations associated with the Schrodingcr equation when using a tensor product of a discrete Fourier ex- 
pansion in space, and a linear interpolation in time. 
The discrete variational equations are given by 



ft 2 fc 2 At 
24tt 2 



At 
24^ 



N/2+j 



forj = -N/2 + !,...,-!, 
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2tt 
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for j = 0,...,N/2-l, 



for j = 0,...,JV/2-l, 



This system of (27V + 3)-equations, allow us to solve for ^fe +1 }^f^Ar/2 an d A i from initial data, 

{v 1 ^ 1 jV l k~ 1 }k= 2 -N/2 anc ^ i^fe' ^fc}fc-Af/2- sucn > this system of equations are an example of a spectral in 
space, second-order in time, pseudo spectral variational integrator for the time-dependent Schrodinger 
equation. The expressions for the variational integrator for the time-independent Schrodinger equation, 
which has spectral accuracy in space, are given by 

N/2+j 

h 2 k 2 ~v^- E' V n -jii n - Xij , for j = -N/2 + 1,..., -1, 

n=-N/2 
N/2+j 

h 2 k 2 ij = — E Vj-nVn — ^Vj , for j = —N/2 + 1, . . . , —1, 

n=-N/2 
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N/2 

h 2 k 2 i, 3 =- Yl Vn-jin - A% , for j = 0, . . . , N/2 - 1, 

n=j-N/2 
N/2 

tfk 2 ^ = - J2' Vj-nVn - Xvj , for J = 0, ... , N/2 - 1, 

n=j-N/2 
N/2 



Vn/2 = — / , v n _ N / 2 V n - AV N / 2 ■ 
n=0 



2 

h 2 k 2 



2 



N/2 

VN/2 = — ^ ^N/2-rfin - At) at/2, 



n=0 
JV/2 

2tt 



1 \ - 



k=-N/2 
V-N/2 — VN/2 ! 
V-N/2 — VN/2 ■ 

As mentioned previously, it is possible generalize this approach to construct a fully spectral variational 
integrator in space-time, using Chcbyshcv polynomials to interpolate in time the coefficients of the discrete 
Fourier expansion used in the spatial interpolation. The computational cost of implementing such a scheme 
would be significantly higher, since this would require all the spatio-temporal degrees of freedom to be solved 
for simultaneously. 

9. Conclusions and Future Work 

We have introduced the notion of a generalized Galerkin variational integrator, which is based on the idea 
of appropriately choosing a finite-dimensional approximation of the section of the configuration bundle, and 
approximating the action integral by a numerical quadrature scheme. 

In contrast to standard variational methods, that are typically formulated in terms of interpolatory 
schemes parameterized by values of field variables at nodal and internal points, generalized Galerkin meth- 
ods utilize function spaces that can be generated by arbitrary degrees of freedom. This allows the intro- 
duction of Lie group methods, and their symmetry reduction using discrete Euler-Poincare reduction, as 
well as multiscale, and pseudospectral methods. Nonlinear approximation spaces allow the construction of 
spatio-temporally adaptive methods, which are better able to resolve shocks and other kinds of localized 
discontinuities in the solution. 

It would be interesting to compare the performance of pseudospectral variational integrators with tradi- 
tional pseudospectral schemes to see if any additional benefits arise from constructing pseudospectral schemes 
using a variational approach. More interesting still would be the comparison for fully spectral methods, since 
both variational and non- variational methods would achieve spectral accuracy, and it would make a particu- 
larly compelling case for variational integrators if their advantages persist even when compared to numerical 
methods with spectral accuracy. 

Most mesh adaptive methods use the principle of equipartitioning the error of the numerical scheme over 
the mesh elements to obtain moving mesh equations. These methods rely on a posteriori error estimators 
that are related to the norm in which the accuracy of the numerical method is measured. While adaptive 
variational integrators exhibit an equipartitioning principle, in the sense that the discrete conjugate mo- 
mentum associated with the horizontal variations are preserved from element to element in each connected 
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component of the domain, it would be interesting to carefully explore the question of whether this can be 
understood as arising from error equipartitioning with respect to a geometrically motivated error estimator. 

While we have only discussed the application of multiscale variational integrators to the case of ordinary 
differential equations, it would be natural to consider their generalizations to partial differential equations, 
whereby the multiscale shape functions are obtained through well-resolved solutions of the cell problem, 
as in the case with multiscale finite elements (see, for example, Hou and Wu [1999]). In general, short- 
term simulations at the fine scale can be used to construct appropriate shape functions to obtain generalized 
Galerkin variational integrators at a coarser level, through the use of principal orthogonal decomposition and 
balanced truncation, for example. This is consistent with the coarse-fine computational approach proposed 
in Theodoropoulos et al. [2000] , or the framework of heterogeneous multiscale methods as proposed in E and 
Engquist [2003]. 

A natural generalization would be to consider wavelet based variational integrators, as well as schemes 
based on conforming, hierarchical, adaptive refinement methods (CHARMS) introduced in Grinspun et al. 
[2002] and further developed in Krysl et al. [2003]. 
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